function ppoly(a,b,m,n)
% a is the left  endpoint, input.
% b is the right endpoint, input.
% m is the number of sets of curves, input.
% n is the degree of polynomial approximation, input.
h = (b - a)/m/n;
X = a:h:b;
Y = f(X);
Xp(1) = a;
Yp(1) = f(a);
for k = 1:m,
  Va = X(n*(k-1)+1);
  Vb = X(n*(k)+1);
  VM = ceil(200/m);
  Vh = (Vb-Va)/VM;
  Vab= X(n*(k-1)+1:n*(k)+1);
  Vyy= f(Vab);
  Px = polyfit(Vab,Vyy,n);
  Vx = Va+Vh:Vh:Vb;
  Vy = polyval(Px,Vx);
  Xp = [Xp,Vx];
  Yp = [Yp,Vy];
end
h = (b-a)/200;
X1 = a:h:b;
Y1 = f(X1);
Z = zeros(X);
plot(X,Y,'or',X,Z,'+r',X1,Y1,'-g',Xp,Yp,'-r');
hold on;
h0 = (b - a)/m;
X = a:h0:b;
Y = f(X);
for k = 1:m+1,
  X0 = X(k);
  Y0 = f(X(k));
  plot([X0 X0],[0 Y0],'-r');
end
hold off;
